
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Test: Exogenous Bidders Participation       %
%               Monte Carlo Experiment              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

L2=100; % number of auctions with N=2
L3=100; % number of auctions with N=3
L=L2+L3;
N=[2,3];
alpha=0.50:0.10:0.80; % quantile levels alpha
psi2=ones(L2,1)*(N(1)*(alpha.^(N(1)-1))-(N(1)-1)*(alpha.^N(1))); % Psi function given N=2
psi3=ones(L3,1)*(N(2)*(alpha.^(N(2)-1))-(N(2)-1)*(alpha.^N(2))); % Psi function given N=3
psi = [psi2;psi3];
m=size(alpha,2);

nsim=10000; % Number of simulations
B=500; % Bootstrap replications
toler=0.1; % Parameter of tolerance for estimation
Rej1=0; % auxiliary sum for rejection probability
Rej5=0; 
Rej10=0; 

% Simulation 

for i=1:nsim
    
       % Data Generating Process
       eps2=normrnd(0,1,L2*N(1),1); % Regression Errors N=2
       eps2=reshape(eps2,N(1),L2);
       eps2=sort(eps2,'descend');
       eps2=eps2(2,:);
       eps3=normrnd(0,1,L3*N(2),1); % Regression Errors N=3
       eps3=reshape(eps3,N(2),L3);
       eps3=sort(eps3,'descend');
       eps3=eps3(2,:);
       
       x2=[ones(L2,1),normrnd(0,1,L2,1),normrnd(0,1,L2,1)]; % Covariates N=2
       x3=[ones(L3,1),normrnd(0,1,L3,1),normrnd(0,1,L3,1)]; % Covariates N=3
       x=[x2;x3];
       p=size(x,2);  
       
       gammai=ones(p,1); % Initial guess for gamma
       
       wb2=x2*gammai+eps2'; % Winning bids N=2
       wb3=x3*gammai+eps3'; % Winning bids N=3
       wb=[wb2;wb3];

       % Computation of the Exogenous Bidder Participation Test Statistic

   for k=1:m
       
       [obj0,res0,~]= MM(psi(:,k), p, wb, x, toler, gammai); % Estimation under H0
       [obj2,res2,~]=MM(psi2(:,k),p,wb2,x2,toler,gammai); % Estimation under H1, sample N=2
       [obj3,res3,~]=MM(psi3(:,k),p,wb3,x3,toler,gammai); % Estimation under H1, sample N=3
       res0k(:,k)=res0; % residuals under H0
       res2k(:,k)=res2; % residuals under H1, sample N=2
       res3k(:,k)=res3; % residuals under H1, sample N=3
       M(k,:)=obj0-(obj2+obj3); 
   end
   
       Mind=sum(M); % Test statistic

       Mstar=zeros(m,B);
   
       % Random Weighting Bootstrap of the Exogenous Bidder Participation Test Statistic  
   
    for b=1:B
        
        w = [mnrnd(L2,ones(L2,1)/L2)';mnrnd(L3,ones(L3,1)/L3)']; % weights
                
        for k=1:m
            res0w=res0k(:,k).*w; % residuals under H0 weighted by w
            res2w=res2k(:,k).*w(1:L2,:); % residuals under H1, N=2, weighted by w
            res3w=res3k(:,k).*w(L2+1:L,:);  % residuals under H1, N=3, weighted by w 
            obj0w=sum(psi(:,k).*res0w)-sum(res0w(res0w<0)); % Objective function under H0 weighted by w
            obj2w=sum(psi2(:,k).*res2w)-sum(res2w(res2w<0)); % Objective function under H1, N=2, weighted by w
            obj3w=sum(psi3(:,k).*res3w)-sum(res3w(res3w<0)); % Objective function under H1, N=3, weighted by w
            Mw(k,:)=obj0w-(obj2w+obj3w);
            % Weighted objective functions
            [obj0wb,~,~]=MMw(psi(:,k),p,wb,x,w,toler,gammai);
            [obj2wb,~,~]=MMw(psi2(:,k),p,wb2,x2,w(1:L2,:),toler,gammai);
            [obj3wb,~,~]=MMw(psi3(:,k),p,wb3,x3,w(L2+1:L,:),toler,gammai);
            Ms(k,:)=obj0wb-(obj2wb+obj3wb);
        end
        
        Mstar(:,b)=sum(Ms-Mw); % bootstraped test statistic
        
    end
    
    % Critical Values
    c99=quantile(Mstar',0.99);
    c95=quantile(Mstar',0.95);
    c90=quantile(Mstar',0.90);
    
    % Computation of rejection probability 
    if (Mind>c99)
    Rej1=Rej1+1;
    end

    if (Mind>c95) 
    Rej5=Rej5+1;
    end

    if (Mind>c90)
    Rej10=Rej10+1;
    end
    
    i
end

% Rejection probability
Rej1=Rej1/nsim
Rej5=Rej5/nsim
Rej10=Rej10/nsim

